% This practical we will work with a single subject's data from a blocked
% finger tapping experiment (CTF data courtesy of Matt Brookes), and
% perform a time-frequency analysis in the beta band using a time-wise GLM
% in source space.   
 
% Data can be downloaded from: 
% www.fmrib.ox.ac.uk/~woolrich/ft_subject1_data.tar.gz

%%%%%%%%%%%%%%%%%%
%% SETUP THE MATLAB PATHS
% make sure that fieldtrip and spm are not in your matlab path

global OSLDIR;
    
tilde='/home/mwoolrich';
osldir=[tilde '/Desktop/osl1.2.beta.11'];

tilde='/Users/woolrich';
osldir=[tilde '/homedir/matlab/osl1.2.beta.17'];

addpath(osldir);
osl_startup(osldir);

%%%%%%%%%%%%%%%%%%
%% INITIALISE GLOBAL SETTINGS FOR THIS ANALYSIS

testdir=[tilde '/Downloads'];
testdir=[tilde '/homedir/matlab/osl_testdata_dir'];

datadir=[testdir '/ctf_fingertap_subject1_data']; % directory where the data is

workingdir=[datadir]; % this is the directory the SPM files will be stored in

cmd = ['mkdir ' workingdir]; unix(cmd); % make dir to put the results in

cd(workingdir);

% Set up the list of subjects and their structural scans for the analysis 
% Currently there is only 1 subject.
clear spm_files;
spm_files={[workingdir '/dsubject1.mat']};
         
structural_files = {[workingdir '/subject1_struct.nii']};      

cleanup_files=0; % flag to indicate that you want to clean up files that are no longer needed

%%%%%%%%%%%%%%%%%%%%
%% load the file in and take a look at the SPM object.

subnum = 1;             
D = spm_eeg_load(spm_files{subnum});

% look at the SPM object. Note that it is continuous data, with 232000 time
% points at 250Hz. We will epoch the data later.
D

%%%%%%%%%%%%%%%%%%%
%% check data in oslview

oslview(D);

%%%%%%%%%%%%%%%%%%%
%% setup osl

oat=[];
oat.source_recon.D_continuous=spm_files;
oat.source_recon.conditions={'Undefined'};
oat.source_recon.freq_range=[13 30]; % frequency range in Hz
oat.source_recon.time_range=[300,32*30];
%oat.source_recon.time_range=[300,12*30];
oat.source_recon.method='beamform';
oat.source_recon.method='mne'
oat.source_recon.gridstep=12; % in mm, using a lower resolution here than you would normally, for computational speed
oat.source_recon.mri=structural_files;
oat.source_recon.dirname=[workingdir '/subj1_results_beta'];
do_hmm=1;
if(do_hmm)
    oat.source_recon.hmm_num_states=-1;
    oat.source_recon.hmm_num_starts=2;
    oat.source_recon.hmm_pca_dim=30;
end;

oat = osl_check_oat(oat);

%%%%%%%%%%%%%%%%%%%%%%%
%% run beamformer

oat.to_do=[1 0 0 0];
oat.source_recon.pca_dim=250;
oat.first_level.name='wholebrain';

oat = osl_run_oat(oat);

%%%%%%%%%%%%%%%
%% Establish regressor for continuous time GLM. 
% This should be setup to correspond to the time window
% over which the source recon is carried out.

D=spm_eeg_load(spm_files{1});

% This should be setup to correspond to the same time window
% as the full time window for D

x=zeros(length(D.time),5);

block_length=30; %s
block_order=[5 5 5 5 5 5 5 5 5 5 4 3 2 1 2 3 1 4 3 4 1 3 2 1 4 4 2 1 3 3 4 1 4 3 1 2 1 2 3 4 3 4 1 2 3 4 1 2];

% [Left, Right, Rest, Both, Rest_at_start]
% [  1     2      4     8     16 ]
% figure;plot(D.time,squeeze(D(1,:,:)))
% emacs ~/vols_data/From_Nottingham_with_Love/JRH_MotorCon_20100429_01_FORMARK.ds/MarkerFile.mrk 

tres=1/(D.fsample);
tim=1;
for tt=1:length(block_order),    
    x(tim:tim+block_length/tres-1,block_order(tt))=1;
    tim=tim+block_length/tres;
end;

figure;plot(D.time,x);

%%%%%%%%%
%% run glm to do regression against known finger tapping regressors

oat.source_recon.dirname=[workingdir '/subj1_results_beta'];
oat=osl_load_oat(oat.source_recon.dirname, 'wholebrain','sub_level','group_level');

oat.first_level.time_moving_av_win_size=0.5;%sec
oat.first_level.cope_type='cope';

% GLM stuff:
oat.first_level.design_matrix=x';
oat.first_level.contrast{1}=[-1 0 1 0 0]'; % rest-left
oat.first_level.contrast{2}=[0 -1 1 0 0]'; % rest-right
oat.first_level.contrast{3}=[0  0 1 -1 0]'; % rest-both
oat.first_level.tf_hilbert_freq_res=diff(oat.first_level.tf_freq_range);
oat.first_level.tf_method='hilbert';
oat.first_level.tf_downsample_factor=10; 
oat.first_level.time_moving_av_win_size=1;

oat.to_do=[0 1 0 0];

oat = osl_run_oat(oat);

% output niis
S2=[];
S2.oat=oat;
S2.stats_fname=oat.first_level.results_fnames{1};
S2.first_level_contrasts=[1:3]; % list of first level contrasts to output
S2.stats_dir=[oat.source_recon.dirname '/' oat.first_level.name '_stats_dir'];
[statsdir,times]=osl_save_nii_stats(S2);    

% make sure you view results using fslview

contrast_num=3;
runcmd(['fslview ' statsdir '/tstat' num2str(contrast_num) '_2mm &']);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Look at ROI time courses
%% calculate an ROI mask 

oat.source_recon.dirname=[workingdir '/subj1_results_beta.oat'];
oat=osl_load_oat(oat.source_recon.dirname, 'wholebrain','sub_level','group_level');
statsdir=[oat.source_recon.dirname '/' oat.first_level.name '_stats_dir'];

% use FSL maths to threshold to create mask
con=3;
thresh=40;
runcmd(['fslmaths ' statsdir '/tstat' num2str(con) '_2mm -thr ' num2str(thresh) ' ' statsdir '/tstat' num2str(con) '_2mm_mask']);

%%
% view the mask in fslview:
runcmd(['fslview ' statsdir '/tstat' num2str(con) '_2mm_mask &']);

%% look at raw timecourse averaged over ROI

oat.source_recon.dirname=[workingdir '/subj1_results_beta.oat'];
oat=osl_load_oat(oat.source_recon.dirname, 'wholebrain','sub_level','group_level');
statsdir=[oat.source_recon.dirname '/' oat.first_level.name '_stats_dir'];

source_recon_results=osl_load_oat_results(oat,oat.source_recon.results_fnames{1});
results = osl_get_recon_timecourses( source_recon_results, [statsdir '/tstat' num2str(con) '_2mm_mask'] );
time_ind=intersect(find(D.time>=oat.source_recon.time_range(1)),find(D.time<=oat.source_recon.time_range(2)));

figure;plot(results.times,normalise(squeeze(mean(results.source_timecourses(1,:),1))));
% compare to design matrix:
ho;plot(D.time(time_ind),x(time_ind,1),'r');
plot(D.time(time_ind),x(time_ind,2),'g');
plot(D.time(time_ind),x(time_ind,4),'k');
legend('data','left','right','both');
plot4paper('time(secs)','beta power');

%% look at Hilbert Envelope timecourse averaged over ROI

oat.source_recon.dirname=[workingdir '/subj1_results_beta.oat'];
oat=osl_load_oat(oat.source_recon.dirname, 'wholebrain','sub_level','group_level');
statsdir=[oat.source_recon.dirname '/' oat.first_level.name '_stats_dir'];

oat.first_level.mask_fname=[statsdir '/tstat' num2str(con) '_2mm_mask'];
oat.first_level.doGLM=0; % does not fit GLM and will output timeseries that would have been input into GLM
oat.first_level.tf_downsample_factor=10; 
oat.first_level.time_moving_av_win_size=4;
oat.first_level.space_average=0;
oat.first_level.name='roi';

oat.to_do=[0 1 0 0];
oat = osl_run_oat(oat);

stats=osl_load_oat_results(oat,oat.first_level.results_fnames{1});
figure;plot(stats.glm_input_times,normalise(squeeze(mean(stats.glm_input_data,1))));
% compare to design matrix:
ho;plot(D.time(time_ind),x(time_ind,1),'r');
plot(D.time(time_ind),x(time_ind,2),'g');
plot(D.time(time_ind),x(time_ind,4),'k');
legend('data','left','right','both');
plot4paper('time(secs)','beta power');


%%%%%%%%%
%% Look at ROI wideband time-frequency spectogram
% we first will need to redo beamformer with a wider band

oat.source_recon.dirname=[workingdir '/subj1_results_beta.oat'];
oat=osl_load_oat(oat.source_recon.dirname, 'wholebrain','sub_level','group_level');
statsdir=[oat.source_recon.dirname '/' oat.first_level.name '_stats_dir'];
con=3;
oat.source_recon.mask_fname=[statsdir '/tstat' num2str(con) '_2mm_mask'];

% need new OAT dirname for new source recon
oat.source_recon.dirname=[workingdir '/subj1_results_wideband.oat'];

oat.source_recon.freq_range=[4 48];
oat.to_do=[1 0 0 0];
oat = osl_run_oat(oat);

%% now do first level for multiple freqs

oat.first_level.mask_fname=[statsdir '/tstat' num2str(con) '_2mm_mask'];
oat.first_level.tf_method='hilbert';
%oat.first_level.tf_morlet_factor=7;
oat.first_level.tf_num_freqs=40;
oat.first_level.tf_freq_range=[4 48];
oat.first_level.tf_hilbert_freq_res=5;
oat.first_level.doGLM=0; % does not fit GLM and will output timeseries that would have been input into GLM
oat.first_level.tf_time_downsample_factor=1; 
oat.first_level.tf_time_moving_av_win_size=10;
oat.first_level.space_average=0;

oat.to_do=[0 1 0 0];
oat = osl_run_oat(oat);
stats=osl_load_oat_results(oat,oat.first_level.results_fnames{1});

figure;imagesc(stats.glm_input_times,stats.glm_input_frequencies,squeeze(mean(stats.glm_input_data,1))');colorbar;
axis xy;
% compare to design matrix:
ho;
x3=x(:,3);
x3(x3==0)=nan;
plot(D.time(time_ind),10*x3(time_ind),'k','LineWidth',5);ho;
legend('rest');
plot4paper('beta power','time(secs)');

%%%%%%%%%
%% do seed based correlation

oat.source_recon.dirname=[workingdir '/subj1_results_beta.oat'];
oat=osl_load_oat(oat.source_recon.dirname, 'wholebrain');

oat.first_level.name=['seedbased'];
oat.first_level.space_average=0;
oat.first_level.time_moving_av_win_size=0.5;%sec
oat.first_level.cope_type='cope';
oat.first_level.doGLM=1;

% remove unwanted fields
try,oat.first_level=rmfield(oat.first_level, 'design_matrix');catch,end;
try;oat.first_level=rmfield(oat.first_level, 'contrast');catch,end;
try;oat.first_level=rmfield(oat.first_level, 'mask_fname');catch,end;

% seed info
oat.first_level.connectivity_seed_mni_coord=[30 -30 48];
oat.first_level.connectivity_seed_regress_zerolag=1;

oat.to_do=[0 1 0 0];

oat = osl_run_oat(oat);

% output niis
S2=[];
S2.oat=oat;
S2.stats_fname=oat.first_level.results_fnames{1};
S2.first_level_contrasts=[1]; % list of first level contrasts to output
S2.stats_dir=[oat.source_recon.dirname '/' oat.first_level.name '_stats_dir'];
[statsdir,times]=osl_save_nii_stats(S2);    

%% make sure you view results using fslview

contrast_num=1;
runcmd(['fslview ' statsdir '/tstat' num2str(contrast_num) '_2mm']);


